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ABSTRACT 

Relativistic jets originating from Supermassive Black Holes (SBHs) can have a con- 
siderable impact on the Interstcllar/Intcrgalactic Medium (ISM/IGM) within which 
they propagate. Here we study the interaction which a relativistic jet, and the co- 
coon associated with its penetration into the ISM, has on the evolution of a dense 
cloud, placed very near the cocoon's path, by analyzing a series of high-resolution 
numerical simulations, and studying the dependence on jet input power, between 
Pj et = 10 41 — 10 47 erg/sec. The density Probability Distribution Function (PDF) 
within the cocoon can be described in terms of two distinct components, which are also 
spatially distinct: a low- and a high-density component. The former is associated with 
the shocked gas within the internal region of the cocoon, while the latter is associated 
with the outer, shocked region of the cocoon itself. The PDF of the post-shocked re- 
gion is well approximated by a modified lognormal distribution, for all values of Pj e t- 
During the active phase, when the jet is fed by the AGN, the cloud is subject both to 
compression and stripping, which tend to increase its density and diminish its total 
mass. When the jet is switched off (i.e. during the passive phase) the shocked cloud 
cools further and tends to become more filamentary, under the action of a back-flow 
which develops within the cocoon. 

We study the evolution of the star formation rate within the cloud, assuming this is 
determined by a Schmidt-Kennicutt law, and we analyze the different physical factors 
which have an impact on the star formation rate. We show that, although the star 
formation rate can occasionally increase, on time scales of the order of 10 5 — 10 6 yrs, 
the star formation rate will be inhibited and the cloud fragments. The cooling time of 
the environment within which the cloud is embedded is however very long: thus, star 
formation from the fragmented cloud remains strongly inhibited. 

Key words: Galaxies: active - intergalactic Medium - large-scale structure of the 
Universe - methods: numerical. 



1 INTRODUCTION 

One of the most intriguing research areas in contemporary 
extragalactic astrophysics involves the study of the inter- 
play between nuclear Black Holes (BHs) in galaxies, the 
(relativistic) jets which they can produce, and the Interstel- 
lar/Intergalactic Medium (ISM/IGM) within which they 
propagate. Observation and modeling of the propagation of 

jets within the ISM is an impo rtant par t of thi s effort. 

Since the seminal wor ks of IScheuerl l| 19741) and iFalld 
(1991), much effort has been dedicated to the study of the 
propagation of a jet into the Interstellar Medium (hereafter 
ISM), and to its consequences for the detectability of the 
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jet. Relatively less attention has been paid to the impact of 
the j et and the cocoon gen erated by it on an inhomogeneous 
ISM. ISaxton et al.l (|2005l ) have performed numerical simu- 
lations of the interaction of a cocoon with a set of clouds 
embedded within a diffuse ISM, mostly paying attention 
to the evolution of t he jet 's morphology. More recently, 
iKrause fc Alexander! (|2007l . hereafter KA07) have also 
studied the turbulence induced by the interaction of jets 
with cold clouds embedded in an Interstellar/Intergalactic 
medium. In their simulations, they show that the jet is 
Kelvin-Helmholtz unstable, and can shear a cold cloud 
embedded in the ISM/IGM, thus inducing turbulence. 
Their simulation setup is however different from ours, 
mainly because they focus their attention on a small 
portion of the jet to study in detail the interaction with 
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IGM clouds and the resulting turbulence. However, as we 
see in our simulations, turbulence in the cocoon arises 
naturally due to the general dynam i cal expa nsion of the 
cocoon itself. ISutherland fc Bicknelll (|2007d lal) have taken 
into account the presence of large-scale density gradients in 
the ISM distribution, and the ability of a jet to emerge out 
of a galactic disc. They notice that the jet can eventually 
percolate through the inhomogeneous ISM, and emerges 
with a very disturbed morphology. 

All the papers quoted made use of fixed-mesh numerical 
codes: in this paper, we use an Adaptive Mesh Refinement 
(AMR) code to follow in detail the evolution of turbulence 
within the cocoon produced by the jet during its propaga- 
tion within the ISM/IGM. Allowing for a high refinement 
level, we can resolve small turbulent eddies, and study their 
statistical properties, and how they affect the evolution of a 
cloud embedded within the cocoon. We focus our attention 
on star formation within the cloud, and on the evolution of 
the density field within the cocoon, as our main aim is to 
understand how, within the turbulent cocoon produced by 
typical AGN jets at moderately high redshifts (z~ 1 — 2), 
this turbulence can modify the ongoing star formation 
within the cloud. 

The paper is organized as follows. In section [5] we briefly 
describe the numerical and computational set-up. In sec- 
tion [2l we present the numerical techniques and physical 
ingredients to simulate the jet and the ISM. We continue 
in section [3] by presenting the initial configuration and the 
physical parameters space spanned by the simulation, and 
we also discuss issues of numerical resolution. In section [5] 
we describe the results of the simulations and the evolution 
of the cocoon. We consider the evolution during both the 
active phase, while the jet is present, and the subsequent 
passive phase. Then, in section [6] we discuss the density 
probability distribution function (PDF) of the matter 
within the cocoon, and we find that it is well described by 
a modified lognormal function. In section [7] we then study 
the evolution of the Star Formation rate (SFR) within the 
cloud. In section [5] we discuss the implications of these 
simulations, comparing our results with those of similar 
papers. Conclusions are presented in the final section. 
In the following, we adopt cgs units, unless otherwise 
explicitly stated, but we adopt kpc for lengths. The under- 
lying cosmological parameters a re taken fro m the 3-year 
WMAP best fit ACDM model (|Beanl I2006T ): H = 74±^ 
Km/sec/Mpc, Q, m = 0.234 ±0.035, n b h 2 = 0.0223 ± 0.0008. 
The unit of time is taken to be the Hubble time for this 
model, i.e.: to = 1.35 x 10 10 years. G, kb and m p denote 
the gravitational constant, the Boltzmann constant and the 
proton mass, respectivelyp 



2 COMPUTATIONAL FRAMEWORK 

To perform the simulations describ ed here, we have used 
FLASH v. 2.5 (|Frvxell et alj l200d ). a parallel, Adaptive 
Mesh Refinement code which implements a second order, 
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shock-capturing PPM solver. FLASH'S modular structure 
allows the inclusion of physical effects like external heating, 
radiative cooling and thermal conduction (among others). 
The user is permitted considerable freedom in the specifica- 
tion of the refinement criteria, which can be customized to 
reach very high spatial and temporal resolutions in selected 
regions. 

The main purpose of this work is to study in detail the inter- 
action of a relativistic jet and the cocoon which it generates 
with pre-existing clouds, and how it can affect star formation 
within the cloud itself. Ideally, a full 3D simulation should 
have a spatial resolution high enough to resolve both the tur- 
bulent motions within the cocoon and the thermodynamic 
structure of the cloud, until the end of the simulation. The 
former task is important when one realizes that, in addition 
to the direct interaction with the jet, the cloud is also sig- 
nificantly affected by the random interactions with the tur- 
bulent eddies present within the cocoon. The computational 
requirements imposed by this task, however, are prohibitive, 
when one considers that the smallest turbulent cells to be 
resolved should have a size comparable to that of the cloud 
(10 h~ x pc here), and the computational box is between 2 
and 4 x 10 2 times larger. For this reason, we have restricted 
ourselves to 2D simulations, where this resolution can easily 
be reached. 

We have included radiative cooling, described by a 
standard cooling function with half solar metallicity 
^Sutherland fc Dopitall993h . extended to high temperatures 
(T > 10 7 K, see appendix [Bj) . Gravity is also included, while 
thermal conduction and magnetic fields are not. The former 
could possibly be relevant for the evolution of the cocoon, 
although on timescales lon ger than those considered here 
l|Krause fe Alexanderll2007h . Regarding magnetic fields, the 
evidence is that, if present within the diffuse IGM, their 
magnitude is not larger than a few microgauss, thus making 
the IGM a high-/3 plasma, and the magnetic field would then 
not significantly affect the global dynamical evolution. 
In the simulations the jet is modeled as a one-component 
fluid, with a density p jet which is a fixed ratio of the en- 
vironmental density p mv . In order to suppress the growth 
of numerical instabilities at the jet/IGM injection interface, 
we adopt a steep, but continuous and differentiable trans- 
verse velocity and density p rofile previously adopted in sim- 
ulations of jet propagation l|Bodo et al.lll994l : |Perucho et al.l 
12004 l2005h : 



V, 



cosii{(y-yj) a '} 



cosh{(y - ViT 3 } 



(1) 



(2) 



where: ctj = 10 is an exponent which determines the steep- 
ness of the injection profile and rij,n env denote the jet and 
environment number densities. This initial profile is highly 
sheared, and peaks around y.j, with most of the thrust 
n jetVj et concentrated around the center of the profile. The 
presence of a highly sheared injection profile forces the code 
to refine the grid structure, populating the injection region 
with subgrid meshes, thus preventing the formation of nu- 
merical contact instabilities. In all our runs, the injection 
point of the jet is chosen to lie at the midpoint of the left 
boundary. 
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cloud is larger. We have chosen to study the effect of the 
jet on a cloud located very near to its path, because this 
allows us to check to what extent star formation is affected 
by the jet under the most extreme conditions. In a forth- 
coming paper we will study a more realistic setup, where we 
distribute a set of clouds with a realistic mass spectrum, and 
study how star formation changes according to the relative 
position within the cocoon associated with the jet. 
Most of the runs were evolved up to « 10 7 yrs., while the 
jet was active and supplying energy to the cocoon. Only in 
one run (H4) the jet was switched off after I0 r yrs., and the 
further evolution of the system was followed until the cloud 
was completely destroyed. 



Figure 1. Initial configuration for the runs. The plot is a mag- 
nification of the actual simulation box, showing a small region of 
the larger simulation box. The small, dense cloud lying along the 
path of the jet, placed at (0.6, 1.92) (in kpch - 1 ). 

3 SIMULATIONS 

Our simulations are characterized by seven parameters. 
Two of these characterize the diffuse ISM: two 
more the cloud: r c i,M c i, and finally three characterize the 
jet: rij, t)j, Vj. As the SFR is mostly determined by the 
strength of the sh ock within the clo ud, as previous models 
seem to suggest (|Klein et al.l 1 19941 ). we have decided to 
mostly vary the jet's input power Pj e t = O.SAjmiinjVj* 
which, together with the density contrast, determines the 
timescale of the cloud's disruption. We keep most of the 
other parameters fixed: ISM density and temperature at 
1 e~cm~ 3 and 10 7 K, respectively, typical of the hot ISM 
in the central parts of a massive elliptical at high redshift 
[z w 1). Only in one run we have increased the size of 
the box, in order to check that the main features of the 
evolution do not depend on boundary conditions. All runs 
except one have been performed on a relatively small box 
(4 kpc), while in run H2 we have used a box twice as 
large. 

In Table 1 we summarize the main parameters of the 
different runs. The ratio between jet and environment 
density, nj/n erl v is fixed to 2 x 10 -2 in runs using small 
boxes, and decreased to a slightly lower value (1CP 2 ) in 
run H2. The injection region has a width of 100 h^ 1 pc 
in runs in the small box, and twice as much in run H2. 
The parameter Vj in eq. [T] is computed once Pj,dj and 
Pjet are assigned. Finally, we have chosen open boundary 
conditions, so the gas is free to flow out of the simulation 
box. The implication of this is that gas is not allowed to 
re-enter the box, so circulation motions on scales larger 
than the simulation box cannot be reproduced. 



3.1 Initial configuration 

We place a small, dense cloud at a position slightly offset 
from the jet's propagation direction (see Fig. [T] and Ta- 
ble |3H|. In order to keep the same spatial resolution, the 
radius of the cloud is doubled in run H2, where the simula- 
tion box is larger: consequently also the initial mass of the 



3.1.1 Model of embedded cloud 

We have chosen a model for the structure of the embedded 
clouds suggested by the numerical simulations performed by 
iBaek et all (|2005l ). because the physical ingredients of these 
simulations are likely to be representative of the physical 
processes present in the ISM/IGM, These simulations have 
been devised to provide a reasonable model for clouds em- 
bedded within the IGM. Different sources of heating (UV- 
background, the wind and radiative flux from the central 
QSO itself) provide a significant energy input which can 
promote the forma tion of pressure-co nfined clouds through 
thermal instability. iBaek et all (|2005h have shown that, for 
typical IGM density and temperatures, similar to those con- 
sidered in the present paper, the cooling time is a small 
fraction of the dynamical time, and the ISM is prone to 
the development of small, pressure-confined clouds. These 
are almost isothermal, with temperatures near the lower ex- 
treme of the cooling function (T c i « 10 4 K). In the mass 
range 10 2 ' 5 ^ M c i ^ 10 7 Mq they find a relationship be- 
tween rt and the total mass M c i : 

n = AMf i4 (3) 

where we have defined: M ctA = M ci /1O 4 M . If rt is mea- 
sured in pc, we obtain: A = 28.87, /3 = 0.28T0.04. The upper 
limit of this mass range corresponds to the Jeans mass for 
these clouds, implying that they are not self- gravitating. 
As we show in Appendix |XJ a reasonable model for 
these clouds is th e Truncated Isothermal Sphere (TIS 
IShapiro et al.| [T999). Given the cloud's mass M c i , the final 
parameters of the configuration will depend on the back- 
ground ISM thermodynamic state, and we assume that the 
latter is described by an unperturbed ideal equation of state 
with density and temperature p c i and T c i, respectively, as 
appropriate to a high temper ature, low density, fully ion- 
ized plasma (e.g. iPriestl Il987l ) . We denote the solution of 
the TIS equation by F(f), so that the cloud density can be 
written as: p c i = poF(Q , po being the normalization factor 
and £ = r/ro a normalized distance. The cloud's structure 
is entirely specified by assigning ro,po and the truncation 
distance r t . Once the l atter is determ i ned b y the radius- 
mass relation found by S hapiro et all (|l999l ). the two re- 
maining factors are determined by imposing pressure equi- 
librium with the IGM, as shown in Appendix lAl 
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Table 1. Parameters of the simulation runs. Columns from left to right are 
as follows: run label, simulation box size, background density (runs Hl-3) or 
central density (runs NFW), jet/background density ratio, jet power. 



run 


L box (kpch 1 ) 


n env (cm 3 ) 




Wj e t (erg/sec) 


HO 


4 


1 


0.02 


4 x 10 40 


HI 


4 


1 


0.02 


8.61 x 10 41 


H2 


8 


1 


0.01 


1.34 x 10 44 


H3 


1 


1 


0.02 


10 45 


H4 


4 


1 


0.02 


10 46 



Table 2. Cloud parameters. All distances are expressed in kpch . From 
left to right, columns are as follows: simulation box size, x and y coordinates, 
cloud's radius and mass, the latter in Mq. 



L box (kpc/i 1 ) 


X 


y 


*cl 


M ci 


4 


0.6 


1.92 


10 


1.65 x 10 8 


8 


1.2 


3.84 


20 


1.30 x 10 9 



4 STAR FORMATION 



4.1 Numerical resolution 



Although our maximum spatial resolution could allow us 
to resolve subparsec scales, we cannot follow star forma- 
tion in any detail, because this would require the inclu- 
sion of more physics (for instance a very detailed treat- 
ment of molecular cooling) and temporal resolution a few 
orders of magnitude higher than the one we have adopted. 
Instead, we assume that the cloud is converting gas into 
stars at a rate determined by the Sch midt-Kennicutt law 
ISchmidtJ ll9o3l . |l95Sl : iKennicuttJ 1 19981 ): S = AT, 71 , with: 
A = (2.5 ± 0.17) x 10" 4 M yr' 1 kpc -2 , n = 1.4 ± 0.15. 
At any time, we assume that star formation proceeds only 
within those regions of the simulation volume where the fol- 
lowing criteria are satisfied: 

(i) Mass is larger than the Bonnor-Ebert mass; 

(ii) Temperature is less than a prescribed upper cutoff, 
i.e. we assume that star formation is sharply inhibited in 
regions having T^T c = 1.2xl0 4 K. 



The Bonnor-Ebert mass l|Ebertlll955l : iBonnorl Il956h is de- 
fined as the largest mass which a pressure confined, gravi- 
tating cloud can reach before becoming unstable: 



M, 



1.18- 



23-55% M Q 

Pext 



(4) 



where p ext is the pressure at the surface of the cloud, tem- 
perature is expressed in units of 10 4 K, and we have as- 
sumed that an isothermal equation of state applies, so that 
the sound speed is given by: c s = (fcsT /irip) 1 ^ 2 . The tem- 
perature cutoff is often adopted in star formation models 
to exclude regions where the UV flux would be too high 
to permit significant star formation. Our criteria to select 
the star forming region are ver y similar to those adopted by 
ISchave fc Dalla Vecchial |200Sl ). 

From now on, all characteristic properties of the cloud such 
as its mass or size will always be referred to the star forming 
region as just defined. 



As our main goal is that of studying in detail the evo- 
lution of star formation within the cloud, we need suffi- 
ciently high spatial resolution during each run. The outer 
parts of the cloud are stripped due to interaction with the 
jet, and it is in principle difficult to predict how much the 
cloud's size and mass will be reduced during the evolu- 
tion. For this reason, we have applied a refinement crite- 
rion which will enable us to get high resolution even dur- 
ing the late phases, when the star forming region is con- 
siderably reduced. The spatial resolution in FLASH is de- 
termined by three parameters: the number of blocks in the 
initial decomposition, which in 2D simulations is given by: 
Nbi x x Nbiy, the number of mesh cells within each block, 
nb x x nb y , and the maximum allowed refinement level l r - 
In the present work we have chosen: N b i x = N b i y = 20, 
Tibx = n by = 8, and l r — 6. Thus, the smallest spatially 
resolved scale (down to mesh level) along each direction is 

given by: A x = A y = (L box /N b i x ) / (2^ r ^n bx ) , which, 

for L box = 4/i _1 kpc gives: A x = 0.78125 h~ 1 pc. Note that 
the maximum refinement level around the cloud is reached 
already at the start of the simulation, because of the refine- 
ment criterion adopted. The approximate number of mesh 
cells contained within a cloud of radius R c i, n m> can be esti- 
mated as: n m = int (nR c i2/ AA m ), where AA m = A x A y is 
the area of a single mesh. In our case initially R c i ~ 10/i -1 
pc, so we get: n m w 514. A typical block distribution is 
shown in Figure [2] 

By default, FLASH refines the grid at those points where 
one of the components of the second spatial derivative of 
some user-selected quantities, normalized to the square of 
the spatial gradient, exceeds some pre-established threshold 
value. The quantities we check for refinement are density, 
pressure and temperature. This default criterion is sufficient 
to resolve the highly compressed regions on the scale of the 
clouds we are interested in, so we do not add any additional, 
customized refinement criterion. 
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Figure 2. Computational block distribution for two different 
steps of run HI, taken at t = 2.6 X 10 5 , upper and 5.3248 X 
10 5 , lower (time in yrs). Each little square represents a FLASH 
block, and each of these contains 64 mesh cells (not reproduced in 
the figures). Note that the size of the region is 50 h~ 1 pc, a small 
fraction of Lf, ox , and the blocks represented are at the highest 
refinement level. 

5 EVOLUTION OF THE COCOON 

The injection of energy into the ISM/IGM can engender tur- 
bulence, mostly because the jet is supersonic at and near the 
injection point. Most previous work aimed at describing the 
general structure and evolution of the cocoon, however, has 
paid more attention to the global dynamics of the cocoon. 
Turbulence can have a significant impact on the evolution 
of the embedded clouds, and for this reason we study it in 
more detail in the next sections. 



5.1 Active phase: evolution during jet injection 

Soon after the jet enters the ISM/IGM a cavity forms, and 
the gas which has been swept out piles up into a transition 
layer. In Figure [3] we show an output of one of the runs at a 
rather advanced stage. We can easily distinguish an internal 
low-density, high temperature cocoon and a shocked ambient 
gas region, externally bounded by a tangential discontinuity 
with the outer ISM/IGM. One of the most interesting fea- 
tures we find is that these regions are also dynamically very 
different. The region containing the shocked gas is threaded 
by a series of weak transonic shocks, and it has on aver- 
age an expansion motion, while in the cocoon a large-scale 
circulation parallel and opposite to the main stream of the 
jet develops, originating from gas reflected away from the 
region near the hot spot. This circulation produces shear 
motions which then decay into weak turbulence within the 
cocoon. Pre-existing clouds embedded within the ISM/IGM 
are heavily affected by this turbulence. The typical veloci- 
ties of these shearing motions are large but, due to the very 
high temperatures within the cocoon (T~ 10 9 — 10 11 K) , 
they are only moderately transonic. 



As one can notice by inspecting Figure[4] the pressure within 
the cocoon can reach high values, because the temperatures 
are on average very high. Also the region around the ter- 
minal part of the jet, near the hot-spot, is subjected to 
high pressures, which increase steadily until the cocoon's 
expansion is halted by the ram pressure. We notice that the 
highest values of pressure are attained within the shocked 
ambient gas region, mostly driven by the higher density, up 
to 3 orders of magnitude larger than the average density in 
the cocoon. The Mach number is higher near the jet, par- 
ticularly near the injection point, but in the overall region 
(cocoon and shocked gas layer) it never reaches values larger 
than M m 3.-3.5. The dynamical evolution of the hot spot, 
i.e. the region between the tip of the jet and the terminal 
part of the cocoon, is quite interesting. Here, however, we 
will concentrate on those features more directly related to 
the interaction with clouds, leaving to further work a more 
detailed analysis of the features of interest relevant to the 
modeling of the radio jet. 

5.2 Passive evolution 

All runs were stopped when the jet was still active, except 
for run H4, which was continued for about 5 x 10 6 yrs. after 
the jet was switched off. This time was chosen to be t w 10 7 
yrs., within the range of a typical duty cycle for the AGN. 
The cocoon expands up to a maximum radius determined by 
the jet power and the environmental ram pressure. When the 
duty cycle of the AGN is completed, the jet injected power 
decreases very rapidly, and the cocoon is not anymore fed by 
the jet. The evolution is mostly driven by inertial motions: 
a snapshot is shown in Figure [5] The typical velocities are 
still quite high within the cocoon, but the expansion of the 
shocked ambient gas region has already been slowed down 
by the ISM/IGM ram pressure, and it is now decelerating, 
while slowly expanding. Notice that the cloud is threaded 
by arc-like internal shocks within the cloud (Figure [7] left), 
previously noticed in similar simu lations of shock-cloud in- 
teractions (Na kamura et al ] l2006ft . 

The further evolution of the cloud during the passive phase 
shows some features not emphasized in previous papers. A 
continuous, decelerating flow now takes place within the co- 
coon, primarily driven by the inertial motions, because at 
the typical densities of the cocoon cm ) 

the gravitational pull is less than the inertial force. A zero 
total velocity line now separates the cocoon from the ex- 
ternal ISM. In the latter the decelerating, outward directed 
expansion motion changes to a wind as the plasma leaks out 
of the simulation box and the density decreases. A similar 
motion, but directed in the opposite direction, drives the gas 
past the left vertical boundary, i.e. towards the AGN and the 
galaxy bulge, increasing the density within this region. As a 
consequence of the generalized density decrease, the cooling 
time increases even further, although it was already much 
larger than the free fall time-scale. 

This back-flow within the cocoon has some interesting con- 
sequences for the embedded clouds. During the active phase, 
the cloud was already subject to a significant compression 
from the turbulent motions, and, as a consequence, its den- 
sity had also risen. Its temperature was however typically 
higher than 5 x 10 5 — 10 6 K, a level maintained mostly 
by the turbulent dissipation. However, during the passive 
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Figure 3. Density and velocity distribution at t = 8.1 X 10 yrs. for run H4. The logarithmic density scale ranges from I0~ i to 10'' 
e~ cm -3 . The transition region between the cocoon and the external medium is threaded by a series of transonic shock waves (shocked 
ambient gas layer). The high-density enhancement within the cocoon originates out of the stripped material from a cloud initially set 
near (0.6,1.92) h~ 1 kpc, which has been shocked by the jet. 



phase, the cloud is embedded within a continuous decelerat- 
ing stream, which is now more laminar than turbulent (Fig- 
ure [6]). This flow continues to compress the cloud, which 
cools efficiently: its average density also increases, and its 
temperature decreases, a trend which can be clearly seen on 
the right hand side of figure [7] On a longer time scale the 
cloud eventually is completely stripped, due to KH instabil- 
ities. Our simulations do not have the spatial and temporal 
resolution to further follow the fragmentation of this cloud, 
which will be treated in a separate paper. The overall ac- 
tion of the cloud's compression during the active jet phase, 
and of the cooling and further compression during the pas- 
sive phase, can result into an occasional enhancement of the 
region of the cloud where favorable conditions for star for- 
mation are possible, as can be seen from figure [7] The gas 
stripped away by KH instabilities from the cloud tends to 



form filamentary, high density structures, where tempera- 
tures can reach high values (T ~ 10 4 — 4 x 10 5 K) . Star 
formation within these filaments, due to these high temper- 
atures, is thus inhibited, and the cloud eventually is com- 
pletely destroyed. 



6 PROBABILITY DENSITY DISTRIBUTION 

The cocoon represents for the cloud an environment radi- 
cally different from the hot ISM: on average, temperatures 
are higher and densities lower, and it is also in a turbu- 
lent state, permeated by a series of turbulent eddies, which 
produce a random, intermittent series of stresses on the 
cloud. This scenario is more complex than those envisaged 
in previous models of cloud evolution under external shocks 
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Figure 4. Density and pressure distribution for the same output of Figure[3] We show two pressure contours, corresponding to 5.77xl0 9 
(white) and 2.88xl0 10 (dark grey) Kern" 3 . 



IIKlein et al, | QUI iNakamura et al.| [2006). and a model of 
the turbulence within the cocoon is then a prerequisite to 
developing a model of the evolution of clouds within a turbu- 
lent environment. We will develop the latter in a following 
paper: here we will only consider the density PDF within 
the cocoon. 

Simulations of turbulence in the ISM suggest that, on galac- 
tic scales, a lognormal Probability Distribution Function 
(LNPDF) provides a viable description of the distribution 
of density fluctuation s over a large range of densities typical 
of ga l actic ISM (e.g. iPadoan et alj 119971 : IWada fc NormarJ 
I2007l . l200ll ). The physical conditions of the system consid- 
ered in this work, however, are very different from those 
typical of the galactic ISM, where temperatures are consid- 
erably lower (T~ 10 — 10 2 K) and densities higher, so we 
do not expect a priori to find the same PDF. Moreover 
the presence of the jet, which injects energy and momentum 
into an expanding cocoon, results in a non-stationary back- 



ground. In Figure [5] we show the evolution of the density 
PDF for one of the runs (H4). One notices that the PDF is 
bimodal, with two peaks at densities respectively lower and 
higher than the initial ambient density (n e =1 cm -3 ). These 
two peaks correspond to two spatially separated regions: the 
low density distribution arises predominantly within the co- 
coon, while the high density component is associated with 
the gas in the shocked gas region. The continuous and dot- 
ted curves are two least-square fits of the low-density PDF, 
with two different fitting functions: an exponentially trun- 
cated power law (continuous curve), 



(5) 



P tr = A I x I™ exp 
and a modified lognormal (dotted curve) 

o [ (z-<*») 2 

_Bexp 



— bx 



(6) 
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Figure 5. Density field of run H4 at At = 2.45 X 10 5 yrs after the jet has been shut off (i.e. at a time t = 1.00245 X 10 7 yrs since 
the beginning of the simulation). The cocoon and the tangential discontinuity are now dynamically very different, and their boundary is 
marked by a zero velocity contour. 



where: x = Log(p). The quality of the two fits is compara- 
ble: x 2 /d.o.f = 1.76 (truncated power law) and \ 2 /d.o.f = 
1.32 (modified lognormal). A modified lognormal distribu- 
tion is generally expected in any fluid system whose den- 
sity distribu tion is the result of a s e ries of uncorrelated 
shocks (e.g yazqucz-Scmadcni 1994; Padoan et al. | 1 19971: 



iPassot fc Vazquez-Semadenil 1 19981 : iKritsuk et~al"1 l2007h 
Such a PDF has also been found in simulat ions of turbulence 
and g l obal star formation in galactic disc s (Wada & N ormanl 
l200ll : iTasker fc Brvanl [2OO6I ; IWada fc Norm an 2007), thus 



suggesting that it could also arise in stationary sys- 
tems where turbulence is not only artificially forced 
within the simulation box. Note that in the model of 
IPassot fc Vazquez-Semadenil (1998) this particular form for 
the PDF should be independent of the dimensionality of 
the system, although the parameters characterizing the dis- 
tribution will depend on it. It seems reasonable to assume 
that the main driver of this weak turbulence are multiple 
shocks within the cocoon, forming initially from the evolu- 
tion of Kelvin-Helmholtz instabilities at the jet-ISM inter- 
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Figure 6. Density evolution at two different epochs after the jet has been switched off. The system evolves towards homogeneity and 
equilibrium on a temporal scale mostly determined by the decay of the inertial motions. Time is measured in units of Hubble time to, 
so the plot on the left is at an epoch t = 2.835 X 10 5 yrs., and the one on the right plot for t = 1.107 X 10 6 yrs. 



face. These shocks then propagate within the cocoon, and 
are reflected at the internal interface with the shocked gas 
region. 



7 EVOLUTION OF STAR FORMATION RATE 

The influence of the jet/cocoon on star formation within 
the cloud is the result of the competition between different 
physical factors. The impact of the shock and the increased 
pressure within the cocoon result in a overall compression 
of the cloud, which tends to increase its density and specific 
star formation rate. On the other hand, the increase of tem- 
perature and stripping of the outer regions due to KH and 
other instabilities tend to reduce the mass of the cloud, thus 
contributing to a decrease of the global star formation rate. 
The detailed temporal evolution of the cloud is governed by 
the turbulence within the cocoon, and subsequently by the 
evolution of the back-flow. In Figure [9] we show the evolu- 
tion of the specific and global star formation rates during 
the active phase, when the jet is still feeding energy into the 
cocoon. One notices that in runs H0-H2 the specific SFR oc- 



casionally increases, due to the increase of the average den- 
sity of the star forming region. The global SFR also shows 
similar fluctuations, although of smaller amplitude, in all 
these runs. A similar trend is observed for the evolution of 
the mass of the star forming region, Figure [POl some episodic 
increases are seen, in temporal and sign correspondence with 
the fluctuations of the SFRs. However, stripping tends to di- 
minish the mass of the cloud, so we are forced to conclude 
that these episodic increase of the star-forming region can 
only be attributed to occasional cooling, driven by thermal 
instability. This gas is however soon warmed up again by the 
warm, turbulent environment of the cocoon, and the SFR 
then decreases. 

A similar trend is observed in run H4, during the passive 
phase. We have already seen from Figure [7] that the cloud, 
exposed to the back-flow from the jet, becomes more fil- 
amentary and cold. However, we see from Figure [11] that 
its mass decreases dramatically, because the inner cold re- 
gions become less dense, and consequently both the specific 
and the global SFRs diminish. We can thus conclude that, 
in general, the interaction with the jet/cocoon tends to in- 
hibit in the long term star formation within the cloud. This 
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Figure 7. Evolution of density and temperature in the compressed cloud, during the active phase (left), during the which the jet is feeding 
energy into the cocoon, and some time after it has been switched off (right). The contours superimposed on the density maps trace regions 
of equal temperature, and correspond (from the inner- to the outermost contour) to the following temperatures: T = 10 4 , 10 7 andl0 9 K. 



suppression however is not continuous, but seems rather to 
proceed in an intermittent way, and is interrupted by short 
episodes during which the SFR occasionally increases. This 
trend seems to continue during the passive phase, when the 
cloud is completely stripped and its SFR decreases rapidly. 



8 DISCUSSION 

Some features of the evolution of the clouds described in 
this work have been previously found in ea rlier work. Our 
initial setup is very similar to that used bv lMellema et all 
(2002): the initial temperatures of the ISM and clouds are 
the same, and our clouds are also in pressure equilibrium 
with the diffuse phase. We also find that the stripped cloud 
seems to fragment into filamentary structures, whose density 
and specific SFR seem to occasionally increase. Although 
these structures seem to persist, we have shown that their 
global star formation rate tends to decrease. 
Recent observational evidence suggests the star formation 
history in early-type galaxies can be more complex than 
in the standard, monolithic scheme. There is both evi- 



dence of quenching of star formation l|Kavirai et all 120071 ; 
Schawinski ct al. 2006) a nd of recent star formation episodes 
(ISarzi et all 120071 : iKavirai et al.l 120071 : iThomas et al]|2007l ; 
lYi et al.ll2007t ): both could be induced by AGN activity. It 
is then important to understand the physical mechanism of 
jet-ISM interaction, and how star formation can be affected. 
The simulations we have presented in this work represent an 
attempt to introduce some realistic physical ingredients of 
the ISM structure into a numerical model, and to under- 
stand how these affect the evolution of SFR in the cloud 
There are some features which have not yet been included 
into this model. One is related to the quantitative rele- 
vance of this model. Although t he duration and freq uency 
of AGN jets are not well known (jBlundell et al.ll2002l ). typ- 
ical estimates of the active phase of the jet range between 
10 7 - 10 8 yrs, comparable to those adopted in the simu- 
lations presented in the previous sections. Empirical evi- 
dence from modeling of SMBH feeding in the context of 
hierarchical structure formation sugg ests that the AGN phe- 
nomenon pers i sts for several Gyr dMahmood et al.l 120051 : 
Bro mley et al. 1 12004 lHaehneltJ[200l ). or many galaxy dy- 
namical times. Hence AGNs are certainly capable of provid- 
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Figure 8. Density PDF for run H4. The analysis was performed 
only within the cocoon/shocked ambient gas region. The two his- 
tograms show the PDF at two different epochs: t=1.35 X 10 4 yrs. 
(dotted), and t=9.45x 10 4 yrs. (continuous). The curves show two 
best fits for the latter epoch: am exponentially truncated power 
law continuous and a a modified lognormal (dotted). The vertical 
line at n e = 7 is the m aximum density possible for fully radiative 
shocks l|BouQuet et al.l [2000) . Best fit parameters are given by: 
(A,b,n) = (1.141 X 10 s , 14.486,24.612) (eq.0, and: (B,(x),cr, b) 
= (3.66 X IO" 3 , -1.6, 0.498, 2.005) (see eqs.[5]and EJ. 



ing significant feedback over the time-scales relevant to star 
formation and gas accretion in hierarchical galaxy formation 
models. 

Our simulations are restricted to only one episode of 
AGN activity: a further episode of jet injection would cer- 
tainly propagate into a very hot and low-density ISM, be- 
cause our simulations show that the diffuse ISM/IGM within 
the cocoon reaches very high temperatures (T« 10 8 — 2 x 
10 11 K) and low densities (n e w 1CT 2 — lO^ 1 cm -3 ), depend- 
ing on the initial temperature and on the jet power at the 
time of its maximum expansion. Under these conditions, 
the cooling time exceeds the dynamical time by a factor: 
tcooi/tdyn ~ 6.5 x 10 2 — 3 x 10 5 , thus the heated gas does 
not c ool significantly (as already noted by llnoue fe Sasakil 
2001). Thus, a second jet would propagate into an environ- 
ment where the ram pressure will be considerably higher, 
Pis, n ~ 10 6 - 2 x 10 10 Kcm" 3 (the initial pressure in the ISM 
was Pism = 10 7 fc _1 ), and would probably be quenched. We 
will verify these predictions in subsequent work. 
Despite the very high temperature within the cocoon, as 
soon as the jet is switched off, we observe that the shocked 
clouds which are not evaporated are subject to a back-flow 
originating from the gas within the cocoon. As a conse- 
quence, they become filamentary, and their temperature de- 
creases, but we have seen that the mass of the star-forming 
regions within the cloud also diminishes, so eventually the 
global SFR is reduced. In a subsequent paper, we will present 
results of a more realistic simulation where we consider the 
propagation of a jet into a multicloud system, and we will 
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Figure 9. Evolution of the cloud's star formation. The red ver- 
tical bar labels the time when the jet/cocoon first hits the cloud. 
Error bars are computed from the ±lcr errors of the Schmidt- 
Kennicutt law (Kennicutt 1998). Upper plot: Specific star forma- 
tion rate, lower plot: Total star formation rate over the whole 
cloud. 



then be able to address the problem of AGN feedback in a 
more quantitative, statistical fashion. 

Our simulations differ significantly from those of KA07, 
because the turbulence within the cocoon surrounding 
the propagating jet is generated naturally by the inter- 
action of the jet with the ISM/IGM, while KA07 stud- 
ied the turbulence associated with the disruption of pre- 
existing ISM/IGM clouds due to Kelvin- Helmholtz insta- 
bilities at the jet-IGM interface. There are however some 
features which we have not a ddressed in our simulations. 
ISutherland fc Bicknelll (|2007bl ) have simulated the propaga- 
tion of a jet which first emerges out of a gaseous disc. They 
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Figure 10. Evolution of the mass of the cloud's star-forming 
region. 
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Figure 11. Mass and SFR evolution after the switch off of the 
jet (run H4). Time is now measured in yrs. after the switch off 
epoch (arbitrarily fixed at t=10 7 ). 



note that the disc has a strong effect on the energy bud- 
get of the jet and also on its morphology: during its initial 
propagation the jet percolates through the disc, and it needs 
more time to overcome this initial pressure and emerge out 
of the disc into the ISM. Sutherland & Bicknell do not con- 
sider the compression of pre-existing clouds: the main tar- 
get of their work is a detailed study of the morphology of 
the emerging jet, and a comparison with observations. A 
significant gaseous component, often distributed in disc-like 
structures, seems however to be present in relevant amounts 
(10 9 - 1O 1O M ) also in a fraction of early-type galaxies 



9 CONCLUSION 

The simulations we have performed in this work have as a 
target the propagation of the jet into the inhomogeneous 
ISM of an host galaxy. We have restricted our interest to 
the central region of a model elliptical galaxy (r < 4 h~ l 
kpc), and we have studied the propagation of the jet and 
its interaction with a single cloud. The adoption of an AMR 
code like FLASH has proved crucial for studying, in the 
same simulation, both the large-scale properties of the jet- 
ISM interaction and the interaction with a small cloud. We 
have modeled star formation within a small, dense cloud 
assuming that it follows the Schmidt-Kennicutt law, and 
we have studied how the SFR is modified by the impact 
of the jet /cocoon. In general, we find that the SFR could 
be occasionally enhanced, but in the long run the cloud is 
stripped by KH instabilities and its SFR decreases. After the 
jet has been switched off, a laminar back-flow flow develops 
which continues to compress and strip the cloud, until the 
latter loses a significant fraction of its mass. Due to the very 
long cooling time of the cocoon, the cloud is embedded in 
a very hot, low-density medium, and this eventually results 
in suppression of star formation. We note that our results 
are preliminary being restricted to simulations in 2D for a 
single cloud. In further work we will extend our study to full 
3D simulations and to an inhomogeneous protogalaxy. 
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APPENDIX A: 
CLOUDS 



PRESSURE CONFINED 



T he truncated spherica l equilibrium configurations studied 
bv lShapiro et al.f(ll999l ) provide a convenient model for IGM 
clouds, because they contain two key physical ingredients: 
they are isothermal and confined by the pressure of an exter- 
nal medium. These configurations are described by solutions 
of a dimensionless Lane-Emden equation: 



d 



d (In p) 



dC 



-PC 2 



(Al) 



Here (, = r/ro and p = p/po are defined as dimensionless 
distance and density, respectively. The characteristic scales 
are related by: 

1/2 



ro 



( ksTci 



(A2) 



(47rGp ) 1/2 \4TrGpm p po, 

where we have introduced the definition: a — (kBT c i/m p ), 
T c i being the cloud's tem perature. The TIS models studied 
by IShapiro et all 0999) are exact solutions of the initial 
value problem for eq. [XT] subject to two initial conditions at 
the origin: 



P (C = o) = i, 



dp 
d{ 



lC=o= 







(A3) 



The first condition enforces the initial value for the den- 
sity, while the second one selects those solutions which are 
not singular, i.e. it defines a core. These configurations still 
have infinitely extending density profiles decaying as r~ 2 
for r < ro, but in presence of an external pressu re the 
profiles get truncated at a radius rt l|McCreal ll957T ). This 
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and: 




Figure Al. Power-law approximation to the TIS mass profile. 
The continuous curve is the function M(C)> the dotted one the 
approximation eq. IA6I 



external confining pressure pt is then an additional parame- 
ter, which will determine the dimensionless truncation dis- 
tance (t- Thus, these models are completely characterized by 
three parameters: pt , the total cloud mass M c i and radius rt ■ 
From these we can deduce the characteristic density po and 
scale ro, by assuming that clouds are desc ribed as TIS pro- 
files. From eq. (35) of IShapiro etafl (|l999h M c i is expressed 
in terms of the averaged mass as: 

M t =r^V /2 J^L (A4) 



Gm p J (47rp ) 1/2 



where: £ t = r±/ro and we have defined: 
M(C) = [ C dss 2 p{s 



(A5) 



In the interval ^ £ ^ 9 the following expression provides a 
good approximation to eq. lA5l with the density p(s) specified 
by the TIS solution: 

M(C) = AC t 1//3 (A6) 

where the coefficient j3 = 0.28 ±0.04 has been introduced in 
eq. [3] This fit is accurate to 4%, as can also be seen from 

Fig. EH 

Inserting eq. IA6I into eq. IA4I and making use of the rela- 
tionship between ro and po (eq. IA2() to eliminate po, one 
obtains: 

2\ 3/2 

g) ^A = ^ (AT) 



Finally, one arrives at an expression for ro: 



ro 
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J- ~1 



\Gm P c< - (M) 

For the values of (3 = 0.28 ± 0.04 and using for c and A the 
values given above, one obtains: 

(A9) 



ro 



0.248 T° t 32 pc 



po = 2.48 x 10 18 M Mpc~ 



(A10) 



APPENDIX B: COOLING FUNCTION 

We adopt the coolin g function provided by 
ISutherland fc Dopital (| 19931 ). for a metallicity [Fe/H] = -1, 
which should be valid in the temperature interval 
10 4 < T < 10 8 K. However, the shocked IGM plasma 
in our simulations often reaches temperatures higher than 
10 8 K, so we need to extend the cooling function to higher 
temperatures. 

The cooling function for a fully ionized plasma in the 
trans-relativist i c regi me (5 x 10 7 < T < 5 x 10 9 K, see 
IWolfe fc Melial I2006T ) has been computed by ISvenssonl 
1 19821 ). and includes contributions 
vant processes like e 4 



from different rele- 
annihilations, and the related 



bremsstrahlung rates. Note that, in principle, cooling from 
e + e~ annihilations is more than 2 orders of magnitudes 
higher than e~p and e~e~ bremsstrahlung cooling rates. 
However, in a relativistic thermal plasma the production of 
e + e~ pairs starts to be significant at temperatures larger 
than the threshold for 77 -> e+e", i.e. for T ^ 5.93 x 10 9 K. 
There e xist few fits to Sv e nsson 's relativistic cooling 
function. IStepnev fc Guilbertl (|l983l ) provide a fit with 4 
coefficients, each computed in tabulated form for different 
values of the temperature. For the more limited temperature 
interval relevant to our work, t wo approximate formula s are 
available from the literature. |Park fc Ostriken l200ll . eqs. 
A1-A4) provide the following fitting formulae: 
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where ot is the Thomson cross section, a/ the fine-structure 
constant, m the number density of ions, and 8 e = kT e /M e c 2 . 
The bremsstrahlung cooling rate Xbr(T e ) is given by: 
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= -Be pn(1.1230 e ) + 1.2746] 6 e > 1. 

This fit is claimed to be accurate to 2% over the full 
transrelativistic domain. 

A relatively simpler fit has been provided by 



A relatively simpler ni 
iKrause fc Alexander (|2007l ): 
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= 2.05 x 10~ 27 Vt (1 + 4.4 x 10~ 10 T) 
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(Note that we have adopted units of erg ■ err? ■ sec 1 in both 
eq.s. 1 and 2). We show these fits in fig. lBll This latter figure 
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Figure Bl. The adopted cooling function. The continuous line 
denotes the Sutherland & Dopita 1993 cooling function, while the 
dotted line is the approximation we have adopted for T ^ 10 8,5 
K, from Krause & Alexander 2005. 

should be compared with Fig. 9 o f lSvenssonl(|l982l ). bvensson 
label the different curves according to the thickness param- 
eter r ni = riir^R = 2.436 x 10~ 7 n 4 [cm -3 ] i?[pc], where: r e 
is the electron radius and R is a typical size of the jet. It 
is clear that, for all the cases we consider, one has r n i <C 1, 
so we can neglect the contribution from thermally-produced 
pairs to the cooling function, which is also dominated by 
bremsstrahlung in this temperature range. 
Note however that no equilibrium relativistic plasma con- 
figuration is possible f or 8 e > 1 (the curve labeled as "O" 
m Fig. 9 ofi venssonl (|l982f ). corresponding to 9 max — 25, 
or: T m ax = 1-48 x lO 11 -^"). At this critical maximum tem- 
perature essentially all the plasma is converted to pairs, 
which then can leave the volume and annihilate outside. For 
T ^ T m ax we then assume A — * oo, n e — * 0. The plasma is 
left with ions, at a significantly lower temperature. 



